In this script we create barplots to show the distribution of survival and reproduction values between the scenarios per combination of stochastic event frequency and severity and we show the distribution between the number of supplements and the years of supplementation for scenarios with positive population growth.

Setup

## for non-CRAN packages please keep install instruction
## but commented so it is not run each time, e.g.
# devtools::install_github("EcoDynIZW/template")

## libraries used in this script
## please add ALL LIBRARIES NEEDED HERE
## please remove libraries from the list that are not needed anymore 
## at a later stage
library("here")
library("tidyverse")
library("magrittr")
library("ggplot2")
library("viridisLite")
library("gridExtra")
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Data

df_StevS_lambda_pos <- readRDS(file = here("output", "data-proc", "model-proc", "df_StevS_lambda_pos.rds"))
df_sc_StevS_mean <- readRDS(file = here("output", "data-proc", "model-proc", "df_sc_StevS_mean.rds"))

Barplot of all scenarios with positive population growth

Version 1: Absolute numbers

Count the individuals a) Overall

table(df_sc_StevS_mean$Mortality_Juveniles)
## 
##  0.2  0.3 0.36 
##  160  160  500
# 0.20  0.30  0.36 
#  160   160   500 
table(df_sc_StevS_mean$Mortality_Subadults1)
## 
## 0.08 0.19 0.26 
##  160  160  500
# 0.08 0.19 0.26 
#  160  160  500
table(df_sc_StevS_mean$Mortality_Subadults2)
## 
## 0.14 0.24 0.31 
##  160  160  500
#0 .14 0.24 0.31 
#  160  160  500
table(df_sc_StevS_mean$Mortality_Adults)
## 
## 0.02 0.14 0.22 
##  240  240  340
# 0.02 0.14 0.22 
#  240  240  340   
table(df_sc_StevS_mean$Repro)
## 
## 0.53 0.58 0.66 1.06 1.41 3.97 
##  400   80   80   80   80  100
# 0.53 0.58 0.66 1.06  1.41 2.88
#  400   80   80   80    80  100
  1. Scenarios with positive population growth
table(df_StevS_lambda_pos$Mortality_Juveniles)
## 
##  0.2  0.3 0.36 
##  160  160  394
# 0.20  0.30  0.36 
#  160   160   394    
table(df_StevS_lambda_pos$Mortality_Subadults1)
## 
## 0.08 0.19 0.26 
##  160  160  394
# 0.08 0.19 0.26 
#  160  160  394   
table(df_StevS_lambda_pos$Mortality_Subadults2)
## 
## 0.14 0.24 0.31 
##  160  160  394
#0 .14 0.24 0.31 
#  160  160  394    
table(df_StevS_lambda_pos$Mortality_Adults)
## 
## 0.02 0.14 0.22 
##  240  214  260
# 0.02 0.14 0.22 
#  240  214  260   
table(df_StevS_lambda_pos$Repro)
## 
## 0.53 0.58 0.66 1.06 1.41 3.97 
##  294   80   80   80   80  100
# 0.53 0.58 0.66 1.06  1.41  3.97
#  294   80   80   80    80   100

Create a matrix

StevS_mort_RR_mat_all <- matrix(data = c(
  500,160,160,NA, # Baseline, 10%, 25%
  500,160,160,NA, # Baseline, 10%, 25%
  500,160,160,NA, # Baseline, 10%, 25%
  340,240,240, NA, # Baseline, 10%, 25%
  400,80,80,80,80,100 # Baseline, 10%, 25%, Status quo, All chicks
), byrow = F)
StevS_mort_RR_mat_all
##       [,1]
##  [1,]  500
##  [2,]  160
##  [3,]  160
##  [4,]   NA
##  [5,]  500
##  [6,]  160
##  [7,]  160
##  [8,]   NA
##  [9,]  500
## [10,]  160
## [11,]  160
## [12,]   NA
## [13,]  340
## [14,]  240
## [15,]  240
## [16,]   NA
## [17,]  400
## [18,]   80
## [19,]   80
## [20,]   80
## [21,]   80
## [22,]  100
StevS_mort_RR_mat <- matrix( data = c(
  394,160,160,NA,
  394,160,160,NA,
  394,160,160,NA,
  260,214,240,NA,
  294,80,80,80,80,100 
), byrow = F)
StevS_mort_RR_mat
##       [,1]
##  [1,]  394
##  [2,]  160
##  [3,]  160
##  [4,]   NA
##  [5,]  394
##  [6,]  160
##  [7,]  160
##  [8,]   NA
##  [9,]  394
## [10,]  160
## [11,]  160
## [12,]   NA
## [13,]  260
## [14,]  214
## [15,]  240
## [16,]   NA
## [17,]  294
## [18,]   80
## [19,]   80
## [20,]   80
## [21,]   80
## [22,]  100

Plot

par(mar = c(5,6,1,1))
barplot(StevS_mort_RR_mat_all, beside = T, ylim = c(0,700), space = c(0,-0.3),  las = 1,
        col = "grey", cex.axis = 1.5, cex.lab = 1.5,  border = NA)
barplot(StevS_mort_RR_mat, beside = T, ylim = c(0,700), space = c(0,-0.3),  las = 1,
        col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 1.5, cex.lab = 1.5, ylab = "", border = NA, add = T)
axis(1, at = c(1.2, 4.9,8.6,12.3,16.5), labels = c("","","","",""), cex.axis = 1.5)
axis(1, at = c(1.2, 4.9,8.6,12.3,16.5), labels = c("s1","s2","s3","s4","RR"), cex.axis = 1.5, line = 0.2, lwd = 0)
mtext("Count", side = 2, line = 4, cex = 1.5)
box()

Version 2: Percentage

Calculate the percent of the used s1-s4 and reproduction values in scenarios with positive population growth

StevS_mort_RR_mat_perc <- (StevS_mort_RR_mat/StevS_mort_RR_mat_all) * 100
StevS_mort_RR_mat_perc
##            [,1]
##  [1,]  78.80000
##  [2,] 100.00000
##  [3,] 100.00000
##  [4,]        NA
##  [5,]  78.80000
##  [6,] 100.00000
##  [7,] 100.00000
##  [8,]        NA
##  [9,]  78.80000
## [10,] 100.00000
## [11,] 100.00000
## [12,]        NA
## [13,]  76.47059
## [14,]  89.16667
## [15,] 100.00000
## [16,]        NA
## [17,]  73.50000
## [18,] 100.00000
## [19,] 100.00000
## [20,] 100.00000
## [21,] 100.00000
## [22,] 100.00000

Plot

par(mar = c(4,8,1,1))
barplot(StevS_mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3),  las = 1,
        col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 2.5, cex.lab = 3, ylab = "", border = NA)
axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("","","","",""), cex.axis = 2.5)
axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("s1","s2","s3","s4","RR"), cex.axis = 2.5, line = 0.2, lwd = 0)
mtext("% of used values in scenarios", side = 2, line = 6, cex = 3)
mtext("with positive population growth", side = 2, line = 4, cex = 3)
box()
text(x = -2.8, y = -5, labels = "b)", xpd=NA, cex = 3)

#784*467
# In figure description: Please consider that overall, survival values increased by 25% were used more often

Save plot

# png(filename = here("plots", "04_PVA", "Stoch_eve_Supplements_bar.png"), width = 1200, height = 772)
# 
# par(mar = c(4,8,1,1))
# barplot(StevS_mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3),  las = 1,
#         col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 2.5, cex.lab = 3, ylab = "", border = NA)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("","","","",""), cex.axis = 2.5)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("s1","s2","s3","s4","RR"), cex.axis = 2.5, line = 0.2, lwd = 0)
# mtext("% of used values in scenarios", side = 2, line = 6, cex = 3)
# mtext("with positive population growth", side = 2, line = 4, cex = 3)
# box()
# text(x = -2.8, y = -5, labels = "b)", xpd=NA, cex = 3)
# 
# dev.off()

A combinded figure of the barplots for the Baseline and Improvement scenarios and the stochastic event and supplement scenarios

# Before this you have to run Script 06c until row 128
# png(filename = here("plots", "04_PVA", "combined_bar.png"), width = 1200, height = 772)
# par(mar = c(0.4,4.6,1,0) + 0.1,
#     mfrow= c(2,1),
#     oma = c(3, 4, 0, 0) + 0.1)
#     #mai = c(1, 0.9, 0, 0.1))
# 
# barplot(mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3),  las = 1,
#         col = viridis(4, end = 0.95), cex.axis = 2.5, cex.lab = 3, border = NA,
#         ylab = " ")
# mtext("% of used values in scenarios", side = 2, line = 6.7, cex = 3, adj = 2)
# mtext("with positive population growth", side = 2, line = 4.7, cex = 3, adj = 2)
# box()
# text(x = -0.7, y = 90, labels = "a)", xpd=NA, cex = 3)
# #text(x = -2.8, y = -5, labels = "positive", xpd=NA, cex = 3)
# 
# barplot(StevS_mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3),  las = 1,
#         col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 2.5, cex.lab = 3, ylab = "", border = NA)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("","","","",""), cex.axis = 2.5, outer = F)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("s1","s2","s3","s4","RR"), cex.axis = 2.5, line = 0.2, lwd = 0, outer = F)
# #mtext("% of used values in scenarios", side = 2, line = 6, cex = 3)
# #mtext("with positive population growth", side = 2, line = 4, cex = 3)
# box()
# text(x = -0.7, y = 90, labels = "b)", xpd=NA, cex = 3)
# dev.off()

Stochastic event frequency and severity

df_StevS_summ <- df_StevS_lambda_pos %>% 
  group_by(Stoch_Frequency, Stoch_Severity) %>% 
  summarise(Freq = n())
table(df_sc_StevS_mean$Stoch_Frequency, df_sc_StevS_mean$Stoch_Severity)
##       
##        0.05 0.1 0.15 0.2 0.25
##   0.05   41  41   41  41   41
##   0.1    41  41   41  41   41
##   0.15   41  41   41  41   41
##   0.2    41  41   41  41   41
table(df_StevS_summ$Freq)
## 
## 33 34 36 37 
##  5  1  3 11
sum(df_StevS_summ$Freq)
## [1] 714
df_StevS_summ$Freq_perc <- (df_StevS_summ$Freq/41)*100 # value of the number of scenarios with positive population growth per combination of stochastic event severity and frequency / the overall number of scenarios per combination of stochastic event frequency and severity * 100 for percentage
  
range(df_StevS_summ$Freq_perc) # 80.4878 90.2439, the range is not so wide
## [1] 80.4878 90.2439

Version 1: Absolute numbers

# Basic Plot
ggplot(data = df_StevS_summ, aes(x = Stoch_Severity, 
                                 y = Stoch_Frequency, 
                                 fill = Freq)) + 
  geom_tile() + 
  # Colour and borders
  scale_fill_viridis_c(option = "viridis") + 
  geom_hline(yintercept = c(0.025, 0.075, 0.125, 0.175, 0.225), 
             colour = "white", 
             size = 1) + 
  geom_vline(xintercept = c(0.025, 0.075, 0.125, 0.175, 0.225, 0.275), 
             colour = "white", 
             size = 1) + 
  # Axis breaks
  scale_x_continuous(#name = "Stochastic event severity",  
                     breaks = c(0.05, 0.10, 0.15, 0.20, 0.25), expand = c(0, 0)) + 
  scale_y_continuous(#name = "Stochastic event frequency", 
                     breaks = c(0.05, 0.10, 0.15, 0.20), expand = c(0, 0)) + 
  # Theme
  theme(rect = element_blank(), 
        axis.text = element_text(size = 14), 
        axis.ticks = element_line(size = 1), 
        axis.ticks.length = unit(.25, "cm"), 
        axis.title = element_text(size = 16), 
        legend.title = element_text(size = 16), 
        legend.text = element_text(size = 14)) + 
  labs(x = "Stochastic event severity", y = "Stochastic event frequency", fill = "Frequency")

Version 2: Percentage

# Basic Plot
ggplot(data = df_StevS_summ, aes(x = Stoch_Severity, 
                                 y = Stoch_Frequency, 
                                 fill = Freq_perc)) + 
  geom_tile() + 
  # Colour and borders
  scale_fill_viridis_c(option = "viridis") + 
  geom_hline(yintercept = c(0.025, 0.075, 0.125, 0.175, 0.225), 
             colour = "white", 
             size = 1) + 
  geom_vline(xintercept = c(0.025, 0.075, 0.125, 0.175, 0.225, 0.275), 
             colour = "white", 
             size = 1) + 
  # Axis breaks
  scale_x_continuous(#name = "Stochastic event severity",  
                     breaks = c(0.05, 0.10, 0.15, 0.20, 0.25), expand = c(0, 0)) + 
  scale_y_continuous(#name = "Stochastic event frequency", 
                     breaks = c(0.05, 0.10, 0.15, 0.20), expand = c(0, 0)) + 
  # Theme
  theme(rect = element_blank(), 
        axis.text = element_text(size = 14), 
        axis.ticks = element_line(size = 1), 
        axis.ticks.length = unit(.25, "cm"), 
        axis.title = element_text(size = 16), 
        legend.title = element_text(size = 16), 
        legend.text = element_text(size = 14)) + 
  labs(x = "Stochastic event severity", y = "Stochastic event frequency", fill = "Percent") + 
  geom_text(mapping = aes(y = 0.05, x = 0.05), label = "a)", vjust = 4.8, hjust = 5, size = 6) + 
  coord_cartesian(clip = "off")

Save plot

The plot with percent is saved

# ggsave(filename = here("plots", "04_PVA", "Stoch_event_percent.png"),
#        width = 13, height = 8.2, units = "cm")

Number of supplements & years of supplementation

table(df_sc_StevS_mean$Supplement_Time, df_sc_StevS_mean$N_Supplements)
##    
##       0  15  30
##   4  20 200 200
##   7   0 200 200
df_Supp_summ <- df_StevS_lambda_pos %>% 
  group_by(N_Supplements, Supplement_Time) %>% 
  summarise(Freq = n())
df_Supp_summ
## # A tibble: 5 x 3
## # Groups:   N_Supplements [3]
##   N_Supplements Supplement_Time  Freq
##           <dbl>           <dbl> <int>
## 1             0               4    20
## 2            15               4   171
## 3            15               7   174
## 4            30               4   174
## 5            30               7   175

Remove the combination with 0 supplements

df_Supp_summ <- df_Supp_summ[-1,]
df_Supp_summ
## # A tibble: 4 x 3
## # Groups:   N_Supplements [2]
##   N_Supplements Supplement_Time  Freq
##           <dbl>           <dbl> <int>
## 1            15               4   171
## 2            15               7   174
## 3            30               4   174
## 4            30               7   175
table(df_sc_StevS_mean$N_Supplements, df_sc_StevS_mean$Supplement_Time)
##     
##        4   7
##   0   20   0
##   15 200 200
##   30 200 200
df_Supp_summ$Freq_perc <- (df_Supp_summ$Freq/200) *100

range(df_Supp_summ$Freq_perc) # 85.5 87.5, very small range
## [1] 85.5 87.5
df_Supp_summ$cat <- c(1:4)

cols <- c("N_Supplements", "Supplement_Time",  "cat")
df_Supp_summ[cols] <- lapply(df_Supp_summ[cols], as.factor)
# Plot with percent of scenarios
ggplot(data = df_Supp_summ, aes(x = Supplement_Time, 
                                y = N_Supplements, 
                                fill = Freq_perc)) + 
  geom_tile() + 
  # Colour and borders
  scale_fill_viridis_c(option = "viridis") + 
  geom_hline(yintercept = c(0.5, 1.5, 2.5, 3.5), 
             colour = "white", 
             size = 1, na.rm = T) + 
  geom_vline(xintercept = c(0.5, 1.5, 2.5), 
             colour = "white", 
             size = 1) + 
  # Axis breaks
  scale_x_discrete(#name = "Time of supplementation",  
                     breaks = c(4,7), expand = c(0, 0)) +
  scale_y_discrete(#name = "Number of supplements", 
                   breaks= c(0,15,30), expand = c(0,0)) + 
   # Theme
  theme(rect = element_blank(), 
        axis.text = element_text(size = 14), 
        axis.ticks = element_line(size = 1), 
        axis.ticks.length = unit(.25, "cm"), 
        axis.title = element_text(size = 16), 
        legend.title = element_text(size = 16), 
        legend.text = element_text(size = 14), legend.box.spacing = unit(2.5, units = "cm")) + 
  labs(x = "Year of supplementing", y = "Number of supplements", fill = "Percent")+ 
  geom_text(mapping = aes(y = 1, x = 1), label = "b)", vjust = 6.5, hjust = 5.8, size = 6) + 
  coord_cartesian(clip = "off")

Save plot

# ggsave(filename = here("plots", "04_PVA", "Supplements_percent.png"), 
#        width = 13, height = 8, units = "cm")

Session Info

## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()
## [1] "2022-01-06 20:15:50 CET"
#git2r::repository() ## uncomment if you are using GitHub
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                    LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3     viridisLite_0.3.0 data.table_1.13.2 magrittr_1.5      forcats_0.5.0     stringr_1.4.0     dplyr_1.0.2       purrr_0.3.4       readr_1.4.0      
## [10] tidyr_1.1.2       tibble_3.0.3      ggplot2_3.3.2     tidyverse_1.3.0   here_0.1         
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.2        pkgload_1.1.0     jsonlite_1.7.1    showtext_0.9      modelr_0.1.8      assertthat_0.2.1  showtextdb_3.0    blob_1.2.1        cellranger_1.1.0 
## [10] yaml_2.2.1        remotes_2.2.0     sessioninfo_1.1.1 pillar_1.4.6      backports_1.1.10  glue_1.4.2        digest_0.6.25     rvest_0.3.6       colorspace_1.4-1 
## [19] htmltools_0.5.0   pkgconfig_2.0.3   devtools_2.3.2    broom_0.7.4       haven_2.3.1       bookdown_0.20     sysfonts_0.8.1    scales_1.1.1      processx_3.4.4   
## [28] d6_0.1.0.0        farver_2.0.3      generics_0.0.2    usethis_1.6.3     ellipsis_0.3.1    withr_2.3.0       cli_2.0.2         crayon_1.3.4      readxl_1.3.1     
## [37] memoise_1.1.0     evaluate_0.14     ps_1.4.0          fs_1.5.0          fansi_0.4.1       xml2_1.3.2        pkgbuild_1.1.0    tools_4.0.3       prettyunits_1.1.1
## [46] hms_0.5.3         lifecycle_0.2.0   munsell_0.5.0     reprex_0.3.0      callr_3.5.0       compiler_4.0.3    tinytex_0.26      rlang_0.4.7       grid_4.0.3       
## [55] rstudioapi_0.11   labeling_0.3      rmarkdown_2.4     testthat_2.3.2    gtable_0.3.0      DBI_1.1.0         R6_2.4.1          lubridate_1.7.9   knitr_1.30       
## [64] utf8_1.1.4        rprojroot_1.3-2   desc_1.2.0        stringi_1.5.3     rmdformats_0.3.7  Rcpp_1.0.5        vctrs_0.3.4       dbplyr_1.4.4      tidyselect_1.1.0 
## [73] xfun_0.18